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Theoretical description of two ultracold atoms in finite 3D optical lattices using 

realistic interatomic interaction potentials 
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A theoretical approach is described for an exact numerical treatment of a pair of ultracold atoms 
interacting via a central potential that are trapped in a finite three-dimensional optical lattice. The 
coupling of center-of-mass and relative-motion coordinates is treated using an exact diagonaliza- 
tion (configuration-interaction) approach. The orthorhombic symmetry of an optical lattice with 
three different but orthogonal lattice vectors is explicitly considered as is the Fermionic or Bosonic 
symmetry in the case of indistinguishable particles. 



I. INTRODUCTION 

The physics of ultracold quantum gases attracts a 
lot of interest since the experimental observation of 
Bosc-Einstein condensation in dilute alkali-metal atom 
gases [l|, [2[- Besides the exciting physics of ultracold 
gases by itself, a further important progress was the load- 
ing of the ultracold gas into an optical lattice formed 
with the aid of standing light waves |3|t5| . The optical 
lattice resembles in some sense the periodicity of a crys- 
tal potential [6|-|9|], but is practically free of phonons. In 
contrast to real solids the lattice parameters are, in addi- 
tion, easily tunable by a variation of the laser intensity, 
the relative orientation of the laser beams, or the wave- 
length. In the first case the lattice depth, in the other 
cases the lattice geometry can be manipulated. 

Moreover, in the ultracold regime the effective long- 
range interaction between atoms is usually well charac- 
terized by a single parameter, the scattering length |10| . 
This simplifies the investigation of atomic systems in the 
ultracold regime. The effective interaction can be either 
attractive or repulsive, depending on the type of atoms 
involved. While different kinds of chemical elements, 
their isotopes, or atoms in different electronic or spin 
states cover already a wide range of interaction strengths, 
an almost full tunability of the atom-atom interactions in 
ultracold gases is achieved using magnetic Fcshbach reso- 
nances [ll|, EQ ■ Close to the resonance value of the mag- 
netic field the scattering length diverges and the effec- 
tive interaction varies in a wide range, in principle from 
being infinitely strongly repulsive to infinitely strongly 
attractive. This possibility of active control of the in- 
terparticle interaction makes ultracold atoms in optical 
lattices an ideal tool for, e. g., exploring the properties of 
many-body Hamiltonians dcscribingparticles in periodic 
potentials, like the Hubbard model |3|, [l3[. Examples are 
the experimental studies of a Bosonic Mott insulator [J] , 
a Fermionic band insulator [5[ , or the simulation of anti- 
ferromagnetic spin chains in an optical lattice [14| . 

A further important aspect of ultracold quantum gases 
in optical lattices is the in principle arbitrary filling that 
can be realized, while the filling is strongly constrained in 
usual solid-state systems by charge neutrality. Ultracold 
quantum gases thus allow studies far away from the often 



considered half-filling case. One interesting limit is the 
one of very sparsely populated lattices in which few- body 
quantum dynamics can be studied very accurately [151 ] . 
This initiated recently a number of corresponding the- 
oretical studies [la - llq . These investigations are often 
further motivated by the fact that the direct experimen- 
tal control of many parameters together with the high 
degree of coherence that is not destroyed by phonons led 
to proposals to use ultracold quantum gases in optical lat- 
tices in quantum-information applications like quantum 
simulators or even quantum computers Q, |{j, [lj, Il9l - t2l| . 
For such applications a very precise knowledge about the 
microscopic interactions between quantum particles in an 
optical lattice is a prerequisite. 

At the required level of accuracy a description of the 
atoms as simply being trapped in an array of harmonic 
potentials becomes inappropriate. For example, during 
a controlled quantum-gate operation it is usually neces- 
sary to bring two atoms from different sites into contact 
with each other. Atoms in different electronic states may 
be involved in such an operation, as those states may 
encode the two qubit states (|0) and |1}). However, in 
such a case the center-of-mass (cm.) and relative (rel.) 
motions of two atoms even in the same potential well 
do not separate, even not within the harmonic approx- 
imation [22, |23|. This is due to the fact that different 
hyperfme states are usually accompanied by different po- 
larizabilitcs. Thus the atoms in different states experi- 
ence different trapping frequencies even within the har- 
monic approximation. Evidently, such a non-separability 
of cm. and rel. motions always occurs for heteronuclear 
systems (different atomic species) with different masses 
and polarizabilites. Experimental evidence for the corre- 
sponding breakdown of the harmonic approximation for 
a heteronuclear system was given in [24| and theoreti- 
cally confirmed [2fJ [26| . If the considered atoms are not 
tightly bound in the same potential well and thus if the 
multi-well structure of the optical lattice is important, 
there is evidently no separability of cm. and rel. motion, 
even not for identical atoms. 

Already the theoretical treatment of two atoms in an 
optical lattice is a formidable task, if realistic atom-atom 
interaction needs to be considered. While this interaction 
may often be described by a central force, the interaction 



potential stems from laborious quantum-chemistry calcu- 
lations and is thus only numerically given. The transition 
to cm. and rel. coordinates simplifies the problem dra- 
matically, since the interatomic interaction affects only 
the rel. motion and in the case of an isotropic interac- 
tion even only the radial part of it. However, the above- 
mentioned non-separability demands finally to treat the 
full six-dimensional problem. Furthermore, the matrix 
elements describing the interaction with the trapping po- 
tential become more involved, since they do not separate 
in cm. and rel. coordinates. Alternatively, the problem 
may be solved in (absolute) Cartesian coordinates. In 
this case the potential describing the optical lattice sep- 
arates for both particles and, e.g., in the orthorhombic 
case even for the three Cartesian coordinates. However, 
in this case the particle interaction terms restore the non- 
separability of the six-dimensional problem, except for 
very special cases like the r 2n potentials [27| where r is 
the radial coordinate of rel. motion and n is an integer. 
In view of its universality with respect to the interparti- 
cle interaction the treatment in cm. and rel. coordinates 
is favorable. Furthermore, as will be shown in detail be- 
low, the use of Taylor expansions of the optical-lattice 
potential allows for a reasonable efficiency. 

In this work a numerical approach is presented for the 
theoretical treatment of two particles that interact via a 
central (isotropic) interaction potential and are trapped 
in a finite orthorhombic sin - or cos 2 -type periodic po- 
tential. While the main motivation is the treatment of 
ultracold atoms and molecules in optical lattices, the 
code allows also to treat other particles and was, e. g., re- 
cently also employed in a study of electrons and excitons 
in quantum-dot molecules. The uncoupled Schrodinger 
equations for cm. and rel. motion are solved by an ex- 
pansion of the radial parts in B splines and the angular 
parts in spherical harmonics. The coupling is then con- 
sidered by means of configuration interaction (exact di- 
agonalization). The orthorhombic D2h symmetry is fully 
accounted for. 

The approach was already successfully applied in a sys- 
tematic investigation of the effects of anharmonicity and 
coupling of cm. and rel. motion for two atoms in a single 
well of an optical lattice [26[. Together with the exper- 
imental results in (2J| this allowed the conclusion that 
only the inclusion of the effects of anharmonicity and cou- 
pling (and thus deviations from a simple uncoupled har- 
monic model) lead to agreement with experiment. Fur- 
thermore, considering a triple-well potential the optimal 
Bosc-Hubbard parameter were obtained and the range of 
validity of the Bosc-Hubbard model was explored quanti- 
tatively [28| . Such a study is of importance, as it provides 
a link to many-body physics and large lattices within the 
most popular model for the description of ultracold atoms 
in optical lattices J3J]. It should be emphasized that in 
view of proposed quantum-information applications the 
physics of, e.g., few atoms in double- w ell p otentials is, 
however, already of interest by itself [29l - l34| . The triple- 
well system was on the other hand, e.g., proposed to 



serve as a transistor, where the population of the middle 
well controls the tunneling of particles from the left to 
the right well [H]. 

The paper is organized in the following way. In Sec. |TT] 
the system and its Hamiltonian is introduced. This 
includes the choice of coordinate system in Sec. Ill Al 
the Taylor expansion of the optical-lattice potential in 
Sec III Bt and the consequent expansion of its angular 
part in spherical harmonics in Sec III CI The obtained 
final form of the Hamiltonian and the alternative cos 2 
lattice are discussed in Sees. Ill C 21 and III C 31 respec- 
tively. Sec. IIIII describes the exact diagonalization ap- 
proach with the corresponding Schrodinger equations in 
Sec IIII Al and all matrix elements that have to be calcu- 
lated in Sec IIIIBI The implementation of symmetry into 
the approach is described in Sec lIVI In Sec [V] computa- 
tional details are given. This includes practical aspects of 
the interaction potential in Sec. IV Al basis-set consider- 
ations in Sec. IV Bl an example calculation of rel. motion 
orbitals for a highly anisotropic trap potential including a 
convergence study in Sec. IV CI and practical issues of the 
final exact-diagonalization step in Sec. IV Dl The paper 
closes with a brief summary and outlook in Sec IVII 



II. 



HAMILTONIAN FOR TWO ATOMS IN AN 
OPTICAL LATTICE 

A. System and coordinates 



The Hamiltonian describing two interacting atoms 
with coordinate vectors f\ , ?% that are trapped in a three- 
dimensional optical lattice is given by 



H{n,f 2 ) 
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Uji^+Uin,^) 



with 



n j (r j )=T j (r j ) + V j (r j ) 



(1) 



(2) 



where Tj is the one-particle kinetic energy operator, Vj 
is the trapping potential of the optical lattice for particle 
j, and U is the atom-atom interaction potential. If the 
lattice is formed by three counterpropagating laser fields 
that are orthogonal to each other, the atoms experience 
the periodic potential 



V J (f J )= J2 V^mn^kcCj), 



(3) 



due to the dipole forces, if the laser frequencies are suffi- 
ciently far-detuned from resonant transitions. In Eq. ([3]) 
VI is the potential depth acting on particle j along the 
direction c (= x, y, z) and is equal to the product of the 
laser intensity I c and the polarizability of the particle j. 
Furthermore, k c = 2tt / \ c is the wave vector and A c is the 
wavelength of the laser that creates the lattice potential 
along the coordinate c. 



A direct solution of the Schrodingcr equation with the 
Hamiltonian given in the form of Eq. (JTJ) is complicated, 
since U depends in general on all six coordinates describ- 
ing the two-particle system, even if the atom-atom inter- 
action is central, i. e., U{f\,r%) = W(|fi— f^l) =: u(r) with 
r = |ri —r%\. For realistic interatomic interaction poten- 
tials, there is no separability and this leads to very de- 
manding six-dimensional integrals. Therefore, it is more 
convenient to treat the two-particle problem in the cm. 
and rel. coordinates, R and r respectively, defined as 



r\ 



T2 



R = filfi + H2?2 



(4) 



with the dimcnsionlcss parameters fj,\ = mi/(mi + m<i) 
and H2 = rn2/(mi + 771,2) where rrij is the mass of the jth 
particle. The system of two atoms in a 3D space as well 
as the cm. -rel. coordinate system is sketched in Fig. Q] 
The evident advantage of this choice of coordinates is the 
fact that the interaction potential acts only on the rel. 
coordinate r and thus on three instead of six dimensions. 
If spherical coordinates are adopted, a central interaction 
potential u(r) depends even on the radial coordinate r 
only. 



Z 1 


1 

f 7 


pi 

R 


\f 


.2 









r 2 


— ► 



FIG. 1: Two particles 1 and 2 in the absolute and cm. -rel. Carte- 
sian coordinate systems. 

On the other hand, the formulation of the two-particle 
problem in the cm. and rel. coordinates complicates the 
treatment of the trapping potential, because its original 
separability in the absolute Cartesian coordinates fj is 
lost in the cm. -rel. system. Only within the harmonic 
approximation for the trapping potential and for two 
identical atoms in the same internal state there is com- 
plete separability in cm. -rel. coordinates |36| . If the true 
atom-atom interaction is furthermore replaced by a 5- 
function pscudopotential, the corresponding Schrodingcr 
equation possesses an analytical solution for isotropic and 
some anisotropic harmonic traps [37ll38j. However, even 
within the harmonic approximation the separability is 
lost, if the two atoms experience different trapping po- 
tentials. This is the case, if a heteronuclear system or 
two identical atoms in different electronic states are con- 
sidered. In the general case of a treatment beyond the 
harmonic approximation and, especially, if the atoms are 



spread over more than one potential well, there is evi- 
dently no separability at all. 

In order to keep the flexibility with respect to the inter- 
particle interaction and the advantage of its simple han- 
dling by using spherical cm. and rel. coordinates, the 
optical-lattice potential has to be brought into a form 
that is convenient for its numerical treatment in those 
coordinates. This is done in two steps. First, a Taylor 
expansion of the sinusoidal trapping potential (J3j> is per- 
formed in Cartesian cm. and rel. coordinates fSec lIIBl . 
The result is then transformed into spherical coordinates 
using an expansion in spherical harmonics (Sec. Ill C|) . 



B. Taylor expansion of the optical-lattice potential 

The optical-lattice potential for both particles in 
Eq. ([3]) is given in the Cartesian cm. -rel. coordinates 
R c and f c (c = x, y, z) as 

2 

v(R,f) = j2 E v c 

i—l c—x,y,z 

xsin 2 (fc c [i? c + (-iy-V„ s r c ]) (5) 

where the index rji = i + (— l) l ~ was introduced for 
compactness. Using the standard trigonometric relations 
the optical-lattice potential can be rewritten in the more 
suitable form 



V(R,r) 



2E E V;{l + (-ir*sm(2k c R c ) 

i—l c=x,y,z 



x sin(2fc c r c /u r)i ) — cos(2k c R c ) cos(2fc c r c /i, K )] . (6) 

In order to achieve a maximal separation of the coordi- 
nates R and r the trigonometric functions in Eq. (|6]) are 
expanded in Taylor series around the origin of the cm. 
and rel. coordinates, 
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sm(2k c R c ) sin(2fc c r c/ i,J = E E (2z + l)!(2j + 1)! 

x {2k c r + \2k c ^ 3 f^Rf +l r^ +l , (7) 

(-l) fc+t 
(2fc)!(2t)! 

x (2k c f{2k c ^ 3 f t Rf rf ■ (8) 

In a practical numerical implementation, the infinite 
sum must be truncated. If, for example, the expan- 
sion is restricted up to the (2n)th degree (with the or- 
der n = 1,2,3,...), the infinite summations in Eqs. ([7]) 
and IJHJ) are changed, since the indices fulfill 

2i + l + 2j + 1 < 2n with 

i < n — 1 — j and j < n — 1 , 

2k + 2t< 2n with 

k <n — t and t <n . (9) 



Hence, the optical-lattice potential can be approximated 
by a Taylor expansion of degree (2n) as 
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where the coefficients 
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are introduced for compactness. 

Using Eq. (fTOj) it is possible to split the optical-lattice 
potential according to 

V(R, f) = v c . m . (R) + v re i. (?) + W(R, r) (13) 

into v c m and v re i. that contain all terms depending solely 
on the cm. coordinate and the rel. coordinate, respec- 
tively. The coupling terms between cm. and rel. motion 
are now contained in W. The three components of V are 

1 2 n 

Vc. m .(i?) = -^E E V^^oTcsR? (14) 



= 1 c=x,y,z fc=l 
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Vrel.(r) =-^E 
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(15) 



where we introduced t re i. and t c . m . for the kinetic-energy 
operators of the cm. and the rel. motion, respectively. It 
is worth emphasizing that in the present formulation only 
the truly non-separable terms (represented by products of 
the cm. and rel. coordinates) are left in the coupling term 
W. All separable terms of the optical lattice potential are 
included into the cm. and rel. Hamiltonians h cm . and 
hrei., respectively. 

In fact, there is a specific case in which the optical- 
lattice potential in Eq. (|6]) can be brought into the form 
of Eq. (fT3| and thus the Hamiltonian (fT7| can be ob- 
tained without performing the Taylor expansion. This is 
the case for two identical particles that are both in the 
same state, if they are deposited in a cubic lattice with 
equal intensities and k vectors along each of the spatial 
directions [22(. If all these conditions are satisfied, then 
Eq. ([S]) can be written as 



,(i?) = 2v E sin2 ( fc ^ c ) 



(20) 
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(r) = 2V E s- 2 ^ , (21 
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W(R,f) = -W E sin 2 (fci? c )sin 2 ( ^- ) . (22) 



Noteworthy, the sum of the sin 2 -shaped lattice po- 
tentials for the two individual particles transforms into 
sin -shaped potentials for both the cm. and rel. motion, 
though the one of the cm. motion possesses a different 
periodicity. In fact, this is also true in the here con- 
sidered general case of a hetcronuclear atom pair in an 
orthorhombic lattice. This is easily seen by extending 
the Taylor expansions in Eqs. (fT4|) to (|16() which gives 



W(i?,r) = iE E V ° 
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(16) 
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lim v c . m .(i?) = V V V c s sin 2 (k c R c ), 



(23) 



lim v rcl .(r) = V V V r /sin 2 (fc c r c /x 7;s ), (24) 

s=l c=x,y,z 



Note that the sum ^ that occurs for n = 1 in the second 

fc=i 
term of Eq. (|16p does not indicate an inverse summation 
but its absence, i.e., no sum at all. 

As a result of the Taylor expansion the Hamiltonian 
(fTj) is transformed into the more convenient form 



H(iJ, f ) = h c . m . (R ) + h re i. (?) +W(R,r) (17) 



with 



.(i?)=t c .m.(i?)+V c . n ,(i?), 



(18) 



hrei. (r ) = t rcl . (r ) + Vrel. (?) + u(r) (19) 
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(25) 



Since an analytical solution for the sin 2 -like lattice 
and non-interacting particles exists |39j , those uncoupled 
known solutions of the Schrodinger equations with the 
Hamiltonians in Eqs. (|2"3"|) and (|24l) could be used as a 
basis for solving the coupled problem. However, due to 



the presence of the central interaction potential a trans- 
formation to the spherical coordinate system is advan- 
tageous, since U{r) = U{^/x 2 + y 2 + z 2 ) does normally 
not allow for a simple solution in Cartesian coordinates. 
However, this change of coordinates is inconvenient, since 
the mentioned analytical solutions of the sin 2 potential 
do not split into simple products of the radial and an- 
gular parts. Since IA can in principle have any possible 
functional form, it is more convenient to transform the 
optical-lattice potential into a form that is suitable for a 
calculation in spherical coordinates. This is done in the 
following subsection by an expansion in spherical har- 
monics. 



C. Expansion of the optical-lattice potential in 
spherical harmonics 

1. Auxiliary functions Y^ mt and Y^ mt 

In order to express the optical-lattices potentials v c . m ., 
v rG L, and W in Eqs. (fT3)l . ([15]) , and (|TB)) . respectively, in 
terms of the spherical harmonics Y™, the corresponding 
polynomials in the Cartesian coordinates r* and Fc c have 
to be rewritten as radial part times an angular function, 



i.e., as r* Fjr 



and .R* -F* c (#, </>), respectively. Every 
and 9 and thus F t c can be ex- 



function of the angles 

panded in spherical harmonics according to 
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i?(M) = £ £ YL*iHM) ( 26 ) 
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where the projection coefficients Y^ mt are given as 
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(27) 



the functions F£{0, 4>) that are equal to cos 2 *(</>) sin 2 '(0) 
(c = x), sm 2 *(c/>)sin 2t (0) (c = y), and cos 2 *(6>) (c = z). 
Consider the integral for Yf mt . The derivation of Yf mt is 
simplified by applying the Euler formula for cos ' (<fi) and 
making use of Eq. (1.111) in Ref. [4Q|. The introduction 
of the new integration variable £ = cos (9) transforms 
the integration limits in (f2"Tf from [0, it] of 9 to [—1, 1] for 
£. Clearly, the integral is non-zero only, if the integrand 
is symmetric in the [—1,1] interval. At this point it is 
important to note that the associated Legendre function 
P/ n is even, if the I + \m\ sum is even, and odd other- 
wise. Since the summation index k in Eq. (1.111) from 
Ref. [40( is an integer that lies in the interval < k < 2t, 
the relation — It < m < 2t is valid. Hence, the m index 
is always even. Therefore, the integral is non-zero only 
for even values of I. Additionally, there are, of course, 
the natural restrictions on the I and m coefficients, i.e., 
I > and 1 77i | < I. Another important fact is that the 

functions Pj^f (£) are oscillatory in the interval 

[— 1 , 1] and the symmetry of the integrand causes the con- 
tribution of negative and positive parts to cancel out, 
leading to a vanishing integral. Hence, one more restric- 
tion on Z is I < 2t. Additionally, Eq. (7.132.1) for the inte- 
gral over the associated Legendre function together with 
Eqs. (8.339.2), (8.339.3) and (8.331.1) (all from Ref. |i|) 
were used in the calculations. 

Summarizing all the steps mentioned above and col- 
lecting the indices that do not give trivial zero contribu- 
tions, the analytical form of Yf mt can be given as 
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I, m even, - 2t < m < 2t, \m\ <l, I < 2t . (29) 



In Eq. ([2T| i stands for the imaginary unit, P] n {9) are the 
associated Legendre polynomials and A;_ m is a constant 
prefactor which is defined by 



Ai m — 



1 21 + 1 (I -m)\ 

4tt (l + m)\ 



(28) 



Finally, the integral over the angular arguments is 

■K 2-K 

Jo d£l = J d9 J d(f>sm(6). Due to their different proper- 

o o 
ties, it is useful to distinguish two types of the expansion 

coefficients Yf mt , those with even values of t and those 
with odd t. The latter ones are in the following denoted 
as Yf mt . According to Eqs. (JT4J) and (|15j) only even pow- 
ers of the Cartesian cm. and rel. coordinates occur in the 
Taylor expansions of v c . m . and v re i. and thus only Y^ mt 
has to be evaluated in those cases. However, the coupling 
term W contains additionally odd powers of R c and r c 
and thus requires also the calculation of Y^ mt . 

First, the calculation of the even-order coefficients 
Y< 



The derivation of the Yf mt coefficient is similar to the 
one for Yf mt and results in 



YLt = (-l)**Fmt 



(30) 



with the same constraints for the indices as for Yf mt . 
In order to derive the Yf f coefficients, Eqs. (7.231.1) 
and (8.752.2) from Ref. J4fl| were additionally used. This 
gives 

oi+2/o, _ i\|| l l 2 - x 

Y\ mt = A l0 (-l)^S mfi (l + ^ +1) „ n H + , 



I even, I < 2t , 



(31) 



where 6ij is Kronecker's delta. 

Finally, the odd-order coefficients Yf TOJ , Y| ' ., 
and Yf m A have to be calculated efficiently. They 



contain the Ff 
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Yf mi , and ^f mt will be considered. They contain cos 2j ' +1 (0), respectively. The Y z c mt coefficients can 
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also be calculated analytically. Consider, for example, 
the integral for Yf m ,. Application of the Eulcr for- 
mula for the cos 2 - ,+1 (0) term and use of Eqs. (1.111) 
and (7.231.2) from Ref. leads to 
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The derivation of Y^ m • is also similar to the one of Yf ■ 
and results in 
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with the same constraints on the indices as for Y 
Finally, Yf m ■ is given by 
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n n— £ 



(-irE E q£x i+1 ^' +l 

j=0 1=0 



Yfo^°( 



E *£-i*7 

m=-i,{2} 



YioA (©^)+ E n A/2 r L A/ (e,$) 

A/=-L,{2} 
! 2t I 

Y / T, cc ^s R2kr2t E E Y L^ m (^» 

t=l fe=l i=0,{2}m=-i,{2} 

2fc L \ 

x E E Y rf(6.$) (37) 

L=0,{2} M=-£,{2} / 



2* It 

is found where, e.g., Yl stands for ^ 

Z=0.{2} 1=0,2,4,... 



In 



- Y y 

l Z0j - ^10] 



and Y 



Imj 



for m 7^ is 



Eq. (03 Y^ 
implied. 

As a result of adopting spherical cm. and rel. coordi- 
nates the Hamiltonian 



2. Final form of the Hamiltonian 

The final expression for the optical-lattice potential is 
obtained by inserting F£(Q,(f>) as defined in Eq. (f2"B| for 
the angular part of the polynomials of R and r occurring 
in Eqs. (fig]). ([T5|). and (JT5J). For the three terms 



H(r, 9, <(>, R, 0, $) = h rcl . (r, 0, 0) + h c , n . (iJ, 9, $) 

+ W(r,6,cj),R,@ > $) (38) 



is obtained with 








hc. m .0R,e,a) = 


i 

~2M 


r ^ 2 a i c 2 . m ; 

<9i? 2 i?<9i? i? 2 








+ v c . m .(i?,e,$ 


) (39) 



1 2 n 

.(r,9,$) = --J2 E K s Ec s fi 2t 

s— 1 c—x,y,z k—1 

2k L 

x E E Y c LMfc y L M (9,$), (35) 

L=0,{2} M=-L,{2} 



Viel.(r,0,< 



^E E KTEcsSLr 3 



s — 1 c=x,y,z t—1 

2t I 

x E E Y? mt Y/"(M), (36) 

l=0,{2}m=-l,{2} 



and 



h re l.(r,l 



1 
271 



d 2 2d 
dr 2 r dr 



+ u r 



+ v re i.(r,0,( 



(40) 



In these equations I c . m . and I rc i. are the operators of 
angular momentum, /i = 77117712/(7711+7712) is the reduced 
mass, and M = mi + 7712 is the total mass of the two 
particles. 

The key achievement is that now all terms in the 
Hamiltonian are at most a sum over products of func- 
tions that depend only on a single coordinate, i.e., 
h{R)f2{ < d)h{&)9i{r)g2{0)g 3 (<j>). As a result all required 
integrals are at most products of one-dimensional inte- 
grals. 



I I I I I I I I I 




FIG. 2: (a) The sin 2 (x) function (black solid) together with 
the 2nd- (blue dashes), 3rd- (red dashes), 5th- (blue solid), 
7th- (red solid) and llth-order (green solid) expansion of the 
Taylor series, (b) The cos 2 (a;) function (black) together with 
the 6th- (green), 7th- (red) and 8th-order (blue) expansion of 
the Taylor series. 



3. Finite and cos 2 lattices 

For practical reasons, the infinite Taylor expansion of 
the sin potential has to be truncated in a calculation. In 
such a situation convergence of the results with respect 
to the expansion length is usually aimed for. However, 
in the present case some caution has to be applied in 
such a study and the interpretation of its outcome [26J . 
The problem is due to the fact that the Taylor expansion 
is performed around the origin. Hence the sin 2 lattice is 
best described close to the origin, while the deviations in- 
crease with increasing distance from it. This is illustrated 
in Fig. HJa) in which the sin 2 lattice (along one coordi- 
nate) is compared to Taylor expansions of different order. 
While the 2nd-ordcr expansion works already very well at 
and close to the origin, even the barrier height of the cen- 
tral well is not correctly described. On the other hand, 
the 3rd-order expansion that includes polynomials up to 
6th order and may thus be called sextic potential agrees 
very well with the central well of the sin 2 potential. The 
sextic potential is thus a very good choice for the investi- 
gation of the effects of anharmonicity on (tightly) bound 
states in a single site of an optical lattice [26| . While in- 
creasing the order of the expansion by considering, e.g., 



the 5th- or 7th-order expansion improves the agreement 
with the sin potential further away from the origin, a 
problem occurs. The resulting potential possesses three 
wells, but the outer ones have a depth that differs from 
the correct one, as is also illustrated in Fig. [Ha). As a 
consequence, non-physical resonances may occur due to 
tunneling. Hence, a simple convergence study in which 
the Taylor expansion is expanded order by order is prob- 
lematic. In fact, an even more severe problem is observed 
for all even-order expansions like the already discussed 
2nd-order one. Those expansions tend to — oo for x go- 
ing to either — oo or +oo (or even in both cases as for the 
shown 2nd-order expansion). As a consequence, an infi- 
nite number of unphysical negative-energy states occur. 

In conclusion, the present approach that is based on a 
Taylor expansion of the optical-lattice potential is espe- 
cially suitable for describing finite optical lattices. With a 
judicious choice of the expansion of the sin potential the 
physics of single- triple-, or higher multiple- well poten- 
tials can be well described. For example, the llth-order 
expansion also shown in Fig. HJa) provides a very good 
description of the physics in a triple- well potential [28| . 
Clearly, even an expansion like the 5th-order one may be 
useful, if an asymmetric potential with different depths 
of the wells should be considered. Furthermore, with a 
sufficiently large number of wells even the physics of a 
complete optical lattice can be described in which the 
continuum states (or transitions to them) arc involved. 
In that case convergence with respect to the number of 
wells (and not with respect to the expansion length) has 
to be achieved, since the true continuum is replaced by 
a correspondingly discrctized spectrum. In fact, it may 
also be reminded that in most experiments involving ul- 
tracold atoms in optical lattices there is an additional 
confining potential and thus the whole (relevant) spec- 
trum can be discrete. 

An evident limitation of the sin lattice and its Taylor 
expansion discussed so far is that only finite lattices with 
an odd number of wells can easily be described. Clearly, 
in many situations also finite lattices with an even num- 
ber of potential wells are of interest. The most promi- 
nent example is certainly the double- well potential that is 
also of special interest for quantum-information studies. 
The physics of few atoms in (one-dimensional) doublc- 
well potentials was recently studied, e.g., in [4l| . The 
most straightforward extension of the present approach 
towards such potentials is by considering a Taylor expan- 
sion of the cos 2 (or 7r/2-shifted sin 2 ) potential 

2 

V C °W) = ]T Yl V c s S in 2 (k c c s + ^). (41) 

s— 1 c=x,y,z 

Using trigonometric relations the optical-lattice potential 
can be written in the more suitable form 

-, 2 
V cos (i?,f) = -£ Y, V c s {l + (-iysm(2k c R c ) 

s— 1 c—x,y,z 

x §\n{2k c r c ^r) s ) + cos(2k c R c )cos{2k c r c tir) s )]- (42) 



After derivations similar to the case of the sin potential 
the splitting of the optical lattice into the cm. and rel. 
motion in the Cartesian frame yields 



*r = Y, E v " 

s—l c—x,y,z 

2 n 

v c co I n.(i?) = 2E E v "zZ { 

s — l c—x.y.z k—1 

. 2 n 



'Okcs^c 



»i = rE E ^E C * 



s — l c—x,y,z t—1 



i 2 

W cos (i?,r*) = -E E ^ 



s=l c=x,y,z 



n-l-j 



(-i) 8 E 

n n — £ 



(43) 
(44) 
(45) 
(46) 



^cos r>2k'2t 



E q£x i+1 ^' +1 + e E QM'Vc 



;=o 



t=i fc=i 



where the constant term v§ os appears that was zero for 
the sin 2 case. Equations (|4"3")) - (|4"(3]) are almost analogous 
to Eqs. ([Ii |) -([T6 ]) for the sin 2 -like potential. In Fig.^fb) 
the cos 2 lattice (along one coordinate) is shown together 
with Taylor expansions of different order. In this case, an 
even-order expansion should be used to avoid unphysical 
negative energy states. For example, while the 6th-order 
expansion provides a rather good representation of two 
neighbor wells in an optical lattice, the 7th-order expan- 
sion leads to negative energy states and does in fact even 
not represent the outer potential barriers properly. The 
8th-order expansion is again rather good, but leads al- 
ready to small artificial side minima. 

Since the adopted expansions are independent for the 
three orthogonal directions x, y, and z, the present ap- 
proach is capable of describing two particles in any com- 
bination of different lattices along those three directions. 
While, c. g., in [28[ a true triple well was described using 
the llth-order expansion of sin in one and a harmonic 
(lst-order) expansion in the two other directions, alterna- 
tively an array of 3 x 3 x 3 lattice sites could be described 
equally well. 



III. EXACT DIAGONALIZATION 

A. Schrodinger equations 

After having formulated the optical lattice potential 
in a suitable form, the solution of the eigenvalue problem 
is described in the following. The Schrodinger equation 
with the Hamiltonian of 



H|* i )=£i|* i ) ) (47) 

is solved by expanding *!/ in terms of configurations, 

* i (5,f)=53C? < ,„*„(5 > f). (48) 



The configurations 



(49) 



are products of the eigenfunctions of the Hamiltonians 
of rel. and cm. motions, respectively, i.e., (p and ip are 
solutions of 



and 



hrol. \<fi) 



hem. \4>j) 



rel. 



<Pi) 



1^5 



(50) 



(51) 



Finally, the wavefunctions of rel. and cm. motion that 
we denote as orbitals (in formal analogy to the one- 
particle solutions in electronic-structure calculations) are 
expressed in basis functions <p and ip that are products 
of B splines and spherical harmonics Y™ for describing 
the radial and the angular parts, respectively, 



w(f> = EC'^ 



(52) 



N r N t I 

E E E Z&m r- 1 B a (r)Yr(8, 0) (53) 

a— 1 I — m— — l 



and 



*(«) = E?'^ 



(54) 



N R N l L 

0=IL=OM=-L 



c5^. M i?- 1 B^(ii)y i M (e,$), 



(55) 



In Eqs. (|5"2"1) and (|5"4"1) . wc introduced the compact indices 
a = a, I, m and b = /3, L, M. 

A specific basis set is characterized by the upper lim- 
its of angular momentum Ni and Nl in the spherical- 
harmonic expansions and the numbers N r and Nr of B 
splines used in the expansions in Eqs. (|53j) and (|55[) as 
well as their order fc ro i. (and fc cm .) and knot sequences 
[421 . |43( . The knot sequences define the ranges of r and R 
in which the wave functions arc calculated, the so-called 
box, though it is, in fact, often a sphere as in the present 
case. If the box is chosen sufficiently large for a given 
finite trapping potential, all wavefunctions will have de- 
cayed before reaching the box boundaries. Otherwise, an 
artificial discretization occurs, if a zero-boundary condi- 
tion at the wall of the box is enforced by removing the 
last B spline. In this case an investigation of the conver- 
gence of the results with respect to the box size has to 
be performed. 

The insertion of the expansions for tp in (|53[) and ip 
in (|55|) into the Schrodinger equations flSOJ) and (|5i"j) , re- 
spectively, followed by a multiplication with either ip* or 
ip* (from left) and integration over r or R leads to gen- 
eralized matrix eigenvalue problems of the type 



he, 



t« SCi 



(56) 



Their solutions provide the energies e\ (and e^ m ) as 
well as the coefficients ££*"• (and c^g 1 ')- The latter define 
the rel. and cm. orbitals ip and ip according to Eqs. (|52|) 
and (|54|) . respectively. Generalized eigenvalue equations 
occur due to the non-orthogonality of the B splines. Fur- 
thermore, the explicit consideration of the factors r _1 
and R~ l transforms the radial part of the Schrddinger 
equations into effective one-dimensional ones by remov- 
ing the d/dR and d/dr terms in Eqs. ([39]) and (gDJ). As 
a consequence, the diagonalizations provide in fact the 
solutions rip and Rip from which ip and ip can, of course, 
easily be obtained. Since rip and Rip vanish for r — > and 
R — > 0, respectively, this additional boundary condition 
is implemented by removing the first B spline. Together 
with the corresponding boundary condition at the outer 
box boundaries, the summations in Eqs. ([5!?]) and ([55)1 

change into J2 a =2 an< ^ Sa=2 > respectively. In fact, 
the actual implementation of the code is flexible with re- 
spect to different choices of the boundary conditions at 
the origin of rel. and cm. motions, but in the following 
only the standard use based on the reduced summation 
limits is considered explicitly for reasons of better read- 
ability. 

Once the eigenvectors ip and ip are obtained, a set of 
configurations $ is built according to Eq. (|4U)) . Again, 
insertion of the expansion for \& in Eq. (|48[l into the 
Schrddinger equation (|47p . multiplication with <!>£ (from 
left), and integration over f and R yields the matrix 
eigenvalue equation 



HC; — Si C; 



(57) 



Due to the orthonormality of the orbitals ip and ip also the 
configurations $ are orthonormal. Therefore, the overlap 
matrix is equal to the identity and Eq. (|57p is an ordinary 
eigenvalue problem. 



For convenience, the integrals over B splines and their 
derivatives are denoted as 



*afa3"a' 



= dri 



,^B Q (r)^B Q /(r) 



drf 



dr v 



(62) 



Furthermore, the index A and the orders of the deriva- 
tives fj, (or v) are omitted for /i = (v — 0). For ex- 
ample, one has B aQ / = V>%> a go a i- Additionally, it is re- 
minded that the character a is reserved for rel. motion 
matrix elements and /3 for cm. elements. Hence, the cor- 
responding notation for the cm. integral over B splines 
analogous to Eq. (|62|) is 



^d»pd"(3' 



= dRR 



x d^B p (R)d' J Bp l {R) 



dR}" 



dR» 



(63) 



Since B splines are polynomials, the integrals B can be 
calculated exactly by means of Gauss-Legendre quadra- 
ture. Due to the compactness (finite local support) of the 
B splines the integration limits are in fact finite. If the 
two involved B splines do not possess a common interval 
where both of them are non-zero, the integral vanishes. 
Therefore, only a very limited number of integrals has to 
be calculated and the resulting overlap and Hamiltonian 
matrices are sparse. In the following, all integrals that 
occur in the calculation are discussed individually. 



The overlap matrices between the basis functions ip 
and ip are not equal to the identity matrix, but 



dnvneAW 



»aa w(C '-'mm 



and, similarly, 



Matrix elements 



'b'b' = ^pp'8ll'8mm' 



(64) 
(65) 



In order to set up the matrix eigenvalue problems in 
Eqs. ([55]) and (j5"7| . the corresponding matrices 

t>&> = (£a|hra.|£a'> , *& = (£.|$v) , (58) 

« = (V ; b|hc. n ,|Vv) , *££ = (^b|Vv) (59) 



and 



#«,«' = (*«|H|^> 



(60) 



have to be set up. As already mentioned, the overlap 
matrix elements between configurations are trivial, 

S K , K < = (*«!$«/} 

= (<Pi K ''Pj K \<Pi Kf ''Pj K >) = &„,V S 3 K -3 K > = ^«,«' ■ ( 61 ) 



2. Kinetic energy 

Since the basis functions are a product of a radial 
B spline and a spherical harmonic, the action of the 
kinetic-energy operator on them is straightforwardly cal- 
culated. Using I? ol Y( m (6l» = 1(1 + l)Y( m (0,0) and 
ic. m .^L M (0, *) = L(L + l)r L M (e, $), one finds 

t V a,8L' = _ ^~% 2 aa' Su'S mm ' + —1(1 + 1) B QQ , 5u<5 mm > 
= Y (BfliaSia' +IQ+ l)«aa') *H'*mm' (66) 

for the rel. motion and analogously 

*b',b ; = 2M {^d^d 1 ^ + L(L + l)B^g,J Sll'Smm' 

(67) 
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for the cm. motion. Note, the second equality in Eq. ([66)1 
as well as Eq. (p7|) have to be (slightly) modified, if non- 
zero boundary conditions are applied at the origin and 
the box boundary. 



Making use of Eq. (fTTj). the angular parts for the ma- 
trix elements of the trapping potential can be calculated 
straightforwardly. For the separable (uncoupled) parts 



3. Interparticle interaction 

The matrix elements of the interparticle interaction 
potential are 



SwS„ 



dru{r)B a {r)B a ,{r). 



(68) 



,rcl. 



\H E K S E C «V E 

/t=0,{2} 



s— 1 c— x,y,z i— 1 



x E n mtt K'i 

m t =-l u {2} 



l t I V \ (l t I V 
I m t m -m'J I 



and 



(72) 



The compactness of the B splines turns the semi- 
indefinite integral into a definite one that has to be calcu- 
lated only within a small spatial interval in which B a and 
B a i (and, of course, u(r)) arc simultaneously non-zero. 
In contrast to the case of the B integrals Gauss quadra- 
ture is in this case only exact, if u(r) can be expressed 
in terms of a finite polynomial expansion. In practice, 
the quadrature converges quite well even with a relative 
small number of terms. This is again partly due to the 
fact that it is sufficient, if a polynomial expansion works 
well piccewise, i.e., only within small spatial intervals. 



4- Separable part of the trapping potential 



Using the property Y^ t *{9,(j)) = (-l) mt YT m *((9, <j>) 
the product of two spherical harmonics can be expressed 
as a sum of products between a spherical harmonic and 
the 3j-Wigner symbols, 



Y, 



i*r< 



h,m t 



t 

m t rn m t 
Here, the coefficient 



o)>T' (».*)• (<») 



A^(-l)V (2t+1>(2 t 1XM+1) <™> 

was introduced for com pac tness. The Gaunt coefficient 
may be obtained as [U, |45| 



dnYr(eA)Y^(e 7 4>)Y/-(e^) 



E 






It I It 

m t m m t 







dVLY~ mt ( 



)Y, 



k i V 



i t iv \ mt rn —ml 



k I I' 




(71) 



v^=-\t E k s E ( 



^cos Tm2fc 
s — 1 c=x,y,z k—1 



2k 



E E 



Y 



L k M k k^L k LL> 



L k =0,{2}M k = -L k ,{2} 

L k L V \ (L k L V 



M k M -M' V 



(73) 



are found for the rel. and cm. matrix elements, respec- 
tively. With the aid of Eqs. flU}, dSSJ), tfHH]), and ^ the 
rel. overlap and Hamiltonian matrices in Eqs. (|58[) arc 
obtained. Insertion into Eq. (|56[) and subsequent diago- 
nalization yields the uncoupled cigenenergics and eigen- 
functions of the rel. motion, as discussed above. Analo- 
gously, Eqs. (|55|) . (|57]). and ([73"]) provide the overlap and 
Hamiltonian matrices in Eqs. (|55|) . thus the eigenenergies 
and eigenfunctions of the uncoupled cm. motion can be 
found. 



5. Matrix elements of the coupled Hamiltonian 

Finally, for obtaining the coupled solutions the Hamil- 
tonian matrix elements H K ^i in Eq. (|60p have to be cal- 
culated. Remind, the total Hamiltonian H was written 
as a sum of the uncoupled Hamiltonians of rel. and cm. 
motion, h re i.,h cm ., and the coupling term W (|38[) . Since 
the configurations are build with the eigenfunctions ip 
and ip of the uncoupled Hamiltonians, only the simple 
diagonal contribution 



($ K |h c . m . +h re i.|$*') 

= (<A K V'jJhc.m. +h. Te i.\ipi K ,tpj K ,) 

= (e™ 1 -+e- m -)^,v^,,V ( 74 ) 



is obtained from h ro i. and h c . m .. 

The remaining task is thus the calculation of the ma- 
trix elements that couple rel. and cm. motions, i.e., the 
ones of W. They are given as 
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Y L k M k k A L k LL' f M * M _ M /J f q* J ? ■ (75) 



Despite the fact that Eq. (ff5|) is somewhat lengthy it 
is convenient and practical for computational purposes. 
While in the computer implementation the summations 
are ordered in such a fashion that the numerical efforts 
are minimized, the order given in Eq. (|75[) is more trans- 
parent. 



IV. SYMMETRY OF THE SYSTEM 

The Hamiltonian of two atoms interacting via a cen- 
tral potential and trapped in sin -like or cos 2 -like po- 
tentials that are oriented along three orthogonal direc- 
tions is invariant under the symmetry operations of the 
point group D 2 h- Since the optical- lattice potential is 
chosen along the three Cartesian axes x, y, and z, the 
single particle Hamiltonians in Eq. ^) T-Lj(x,y, z) = 
fij(-x>-y,-z) = Hj{-x,y,z) = Hj(x,-y,z) = 
Hj(x,y,-z) = Hj(-x,-y,z) = Uj(x,-y,-z) = 
Hj(—x,y,—z) arc equivalent. This is a consequence of 
the symmetry elements of the orthorhombic D 2 h group 
that contains besides the identity operation E and the 
inversion symmetry i also three twofold rotations (by the 
angle n) C 2 (x), C 2 (y), and Ci(z) as well as the three 
mirror planes <j(xy), a(xz), and a(yz). The symmetry 
elements are illustrated in Fig. [3] 

The symmetry group D 2 h has eight irreducible repre- 
sentations (see Table IJ): A g , B lg , B 2g , B 3g , A u , B\ u , 
S 2m , B 3u . Clearly, the explicit use of symmetry is ad- 
vantageous, since it reduces the numerical efforts dra- 



C 2 (z) 



O(xz) 




0"Uy) 



' C 2 (y) 



G(yz) 

FIG. 3: The symmetry elements of the two particles interacting 
by a central potential in a sin 2 -like trap. The list is complete with 
the identity clement E added. 



matically as the different irreducible representations can 
be treated independently of each other. This reduces the 
dimensions of the matrices that have to be diagonalizcd 
by approximately a factor of 8 x 8 = 64. Furthermore, 
many integrals vanish for symmetry reasons and have 
thus not to be calculated at all. 

In fact, the D 2 h symmetry is a consequence of the con- 
sidered shape of the potential and thus the symmetry 
of the single-particle Hamiltonians in absolute Cartesian 
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TABLE I: Character table of the D2h point group 



D 2h 


E 


C 2 (z) 


C 2 (y) 


C 2 {x) 


i 


o(xy) 


oixz) 


a(yz) 


A g 


1 


1 


1 


1 


1 


1 


1 


1 


Big 


1 


1 


-1 


-1 


1 


1 


-1 


-1 


B2 9 


1 


-1 


1 


-1 


1 


-1 


1 


-1 


B 3a 


1 


-1 


-1 


1 


1 


-1 


-1 


1 


A u 


1 


1 


1 


1 


-1 


-1 


-1 


-1 


B\u 


1 


1 


-1 


-1 


-1 


-1 


1 


1 


B 2 U 


1 


-1 


1 


-1 


-1 


1 


-1 


1 


Bsu 


1 


-1 


-1 


1 


-1 


1 


1 


-1 



TABLE II: Results of the D2h group operations on the abso- 
lute and spherical coordinates, and the corresponding trans- 
formations of spherical harmonics. Given are the values of 
a, b, c that are multipliers for the x, y, z coordinates, respec- 
tively, and 9' , <j>' that are shifts of the spherical coordinates 
9, <f>, respectively. 



symmetry 


Absolute 


Spherical 


yrn 




(ax, by, cz) 


{9' + 9, cj>' + <j>) 


Y l m {9' + 9,4>' + <t>) 


E 


(1,1,1) 


(0 + 61,0 + 0) 


Yr(9,<p) 


C 2 (z) 


(-1,-1,1) 


(0 + 9,ir + 4>) 


(-i) m YT(o,<i>) 


Ca(y) 


(-1,1,-1) 


(ir - 9,ir - <(>) 


{-l) l+m Y- m (8,4>) 


C 2 (x) 


(1,-1,-1) 


(tt - 0,2tt - (f>) 


{-l) l Y-™{6A) 


i 


(-1,-1,-1) 


(ir-9,-K + <j>) 


{-i) i Yr(8,<t>) 


a(xy) 


(1,1,-1) 


{■k- 9,0 + 0) 


{-i) l+m Yr{9,<t>) 


a(xz) 


(1,-1,1) 


(0 + 0, 2tt - <£) 


(-l) m Y l - m {9A) 


a(yz) 


(-1,1,1) 


(Q + 9,TT-(f>) 


Y- m {6,0) 



coordinates, TLj in Eq. (J2|). Since the atom-atom interac- 
tion u is invariant under all operations in D 2 h, the total 
Hamiltonian H in Eq. ([1]) belongs also to the D 2 h group. 
As may be less transparent on a first glance, but can be 
shown from a complete symmetry analysis, also the rel. 
and cm. Hamiltonians h re i (f) and h c m (R) possess the 
same symmetry as the total Hamiltonian. Therefore, it is 
sufficient to examine the symmetry properties, e.g., for 
the rel. part only. The cm. part has the same properties 
and the ones of the total Hamiltonian can then be de- 
duced from the properties of the direct tensor products. 
In order to use the symmetry when solving the eigen- 
value problems of the uncoupled rel. and cm. motions, 
symmetry-adapted basis functions have to be obtained. 
Since all basis functions adopted in this work are cen- 
tered at the origin and are products of a radial part 
times a spherical harmonic, ([55)1 and (|55|) . the symme- 
try operations affect only the angular part. Therefore, 
linear combinations of the spherical harmonics have to 
be found that transform like the irreducible represen- 
tations of the D 2 h group. The problem of determining 



symmetry-adapted basis functions from the complete set 
of spherical harmonics has a long history, starting with 
the introduction of cubic harmonics for the cubic point 
group [461 ] . Since it appears, however, to be not that 
trivial to find the orthorhombic harmonics in easily ac- 
cessible form, they are given explicitly together with a 
brief derivation. First, the action of the symmetry ele- 
ments of D2h on the spherical harmonics has to be con- 
sidered. The result is shown in Table [XT] that provides 
also the intermediate steps, the result of applying the 
symmetry operations on the Cartesian and the spheri- 
cal coordinates. The most important result is that all 
symmetry operations of D 2 h leave I unchanged, i. e., only 
Y" 1 with identical values of I are transformed into each 
other. As a consequence, the symmetry-adapted basis 
functions are superpositions of spherical harmonics with 
a fixed value of I. 

In view of the eight irreducible representations of D 2 h 
(see Table H| one needs to find the required eight sets of 
orthonormal linear combinations of spherical harmonics. 
In the present case, this is easily achieved using the stan- 
dard projector technique, i. c, by applying the projector 



rE*(6* 



O fc 



(76) 



k=\ 



of the irreducible representation i onto a spherical har- 
monic Y™ . In Eq. (|76]l it is used that all irreducible 
representations in D 2 h are non-degenerate. Furthermore, 
h is the total number of symmetry operations (eight for 
D 2 h), 0^ the operator corresponding to symmetry ele- 
ment k, and Xi(Ofc) the character of symmetry element 
k for the irreducible representation i. While all char- 
acters are listed in Table U the results of the symmetry 
operations on Y{ a are given in the last column of Table [Til 
For example, the application of Ps 3u gives 



p s 3u Y i m = - [1 - (- 1 )'" - (-!)' + (-l)' +m l Y t m 
8 

+ [(-!)' - (-l) ; + m + (-1)™ - l] Yf m . (77) 



Clearly, only the combination of odd values of I and m 
yields in this case a non-zero result and thus a symmetry- 
adapted basis function, 



Pss. vr = 2 ( Y i m - Y r n ) l > m odd • ( 78 ) 



The use of these symmetry-adapted basis functions (su- 
perposition of spherical harmonics instead of a single one 
and restriction on I and m) modifies the wave functions 
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TABLE III: Product table of the D 


2h point group 


<8> 


A g 


Big 


B 2 g 


Ba a 


All 


Biu 


B2u 


Bau 


A a 


A g 


Big 


B 2 g 


Bag 


-'It/ 


Biu 


B2u 


Bau 


Big 


Big 


A g 


Bag 


B2g 


Biu 


Au 


Bau 


B 2 u 


B2g 


B2g 


Bzg 


A a 


Big 


B2u 


Bau 


A u 


Biu 


B Sa 


Bag 


B2g 


Big 


A a 


Bau 


B 2 u 


Biu 


A-u 


A-u 


Jiu 


Biu 


B2u 


Ba u 


A a 


Big 


B 2 g 


Bag 


B\u 


Biu 


- i u 


Ba u 


B2u 


Big 


A g 


Bag 


B 2g 


B2u 


B2u 


Bzu 


A u 


Biu 


B 2g 


Bag 


A B 


Big 


Ba u 


Bz u 


B2u 


Biu 


A v 


Ba g 


B2g 


Big 


Ag 



of the rcl. motion in Eq. (|53|) into 

N r Ni I 

^ = E E E Z&knr- 1 B a (r)9& 

a=l i=0,{2} m=0,{2} 

N T Ni I 

vf 19 = E E E Z&L'- 1 B a {r)9fc 

a=l l=2,{2}m=2,{2} 

N r Ni I 

v! 29 = E E E tfSLr- 1 B*{r)9fa 

a=l l=2,{2}m=l,{2} 

v! 39 - E E E ^ m ^^w+ 

a=l i=2,{2}m=l,{2} 

N r N, I 

<& = E E E ~ct lm r- 1 B a {r)W lm 

a=l i=3,{2}m=2,{2} 

N r Ni I 

^ - E E E tf^r- 1 B a (r)*& 

a=l i=l,{2}m=0,{2} 

JV r N l I 

<#** = E E E tfS^r- 1 B a {r)9+ 

a=l J=l,{2}m=l,{2} 

' = E E E ^ lm r^B a (r)^- 

a=l J=l,{2} m=l,{2} 



•r 
where 



^o + = ^o~=^ (M), 



m = Yr{e^)±Yr m {0^) (m^O) 



79) 
80) 
81) 
82) 
83) 
84) 
85) 
86) 

87) 



is introduced for compactness and a summation index 
/ = i, {2} means I — i, i + 2, . . . Moreover, for the coeffi- 
cients c the superscripts rel. arc omitted for better read- 
ability. Clearly, the consideration of symmetry-adapted 
basis functions for the cm. motion leads to a completely 
analogous modification of Eq. ([55]) . 

Since the D211 point group contains only non- 
degenerate irreducible representations, its product ta- 
ble (showing the result of a direct tensor product be- 
tween a pair of irreducible representations and given in 



Tab. IIIip is straightforwardly obtained and every prod- 
uct corresponds uniquely to one irreducible representa- 
tion. For example, a configuration $ K = ip i 3g ip- 2s trans- 
forms as Big. Clearly, symmetry- adapted configurations 
are straightforwardly constructed from the symmetry- 
adapted rel. and cm. orbitals. 

In the case of indistinguishable atoms, the quantum 
statistics has to be considered. For Fermionic atoms, the 
total wavefunction must change sign under particle ex- 
change, while it must remain the same for Bosons. Parti- 
cle exchange does not influence the cm. coordinate (i. e., 
R — > R or, equivalently, $ — >• $, 9 — > 6), but the rcl. 
coordinate (i. c, r — > -for <fi — >• ir + cf), 6 — >• 7r — 9). Since 
particle exchange corresponds to the symmetry operation 
of inversion (i) for the coordinate, all gerade basis func- 
tions (ifiAg, fB lg , t PB 2g i an d l PB 3g ) are allowed for identi- 
cal Bosons, the ungerade functions (<PA n > fB lu > Pb 2u > an d 
<Pb 3u ) for identical Fermions. The quantum statistics for 
indistinguishable atoms is thus easily taken into account 
and reduces the number of possible orbital combinations 
by factor 2. The straightforward (almost automatic) im- 
plementation of the quantum statistics that leads even 
to a direct reduction of the computational demands is a 
further advantage of the present approach. 



V. COMPUTATIONAL DETAILS 

The theoretical approach presented in this work pro- 
vides an efficient way to treat two interacting atoms in 
an orthorhombic optical lattice. The use of cm. and rel. 
coordinates and the expansion of the basis functions in B 
splines for the radial part and spherical harmonics for the 
angular parts is especially useful for considering realistic 
interatomic (molecular) interaction potentials. Further- 
more, the Bosonic or Fermionic nature of the atoms is 
easily accounted for. However, in the case of strongly 
anisotropic lattice potentials the advantage of the use 
of spherical harmonics that all involved integrals can be 
efficiently and analytically calculated is partly compen- 
sated by their slow convergence. Similarly, the adopted 
exact diagonalization approach for incorporating the cou- 
pling of cm. and rel. motion has the advantage of being 
exact, if converged, but is known to be slowly conver- 
gent. For many experimentally relevant parameters the 
calculations are, therefore, still very demanding and thus 
an efficient implementation of all involved computational 
steps was mandatory. The adequate choice of the basis- 
set parameters adapted to the considered problem can 
improve the efficiency drastically. Thus it is worthwhile 
to at least briefly discuss some technical aspects of both 
the implementation and the choice of basis sets. 



A. Interatomic interaction potential 

In the present approach the interatomic interaction en- 
ters the calculation only in the determination of the rcl. 
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orbitals. In the so far considered case of isotropic inter- 
action potentials the potential influences only the calcu- 
lation of the radial integral (|68)) . Clearly, an extension to 
orientation-dependent interaction potentials (like dipole- 
dipole interaction) is possible by expanding the angular 
part of the interaction potential in spherical harmonics. 
Then the resulting angular integrals can still be solved 
analytically. Since the radial integral (|55|) is calculated 
using quadrature, even a usually only numerically given 
Born-Oppenheimer potential curve can be used in order 
to achieve a realistic description of the interatomic inter- 
action. Clearly, any type of potential can be easily used. 
Only the implementation of the (^-function pscudopoten- 
tial required some care due to its numerically ill-behaved 
nature. 

In the ultracold regime the scattering process is ex- 
tremely sensitive to all details of the complete inter- 
atomic potential curve. In fact, for most experimentally 
relevant alkali atoms, it is impossible to calculate the 
potential curves with sufficient precision. In the zero- 
energy limit the scattering process is fully described by 
the scattering length a sc [10(. Depending on the consid- 
ered atoms, their isotopes, and electronic states a sc can 
have very different values for different systems, ranging 
from a sc 3> (strongly repulsive) via a sc m (almost 
non- interacting) to a sc <C (strongly attractive). An 
important aspect of ultracold atomic gases is the fact 
that many atom pairs possess magnetic Feshbach reso- 
nances, i.e., with a magnetic field the colliding atoms 
can be brought into resonance with a molecular bound 
state. As a consequence, the scattering length becomes 
experimentally tunable [IJ, |47|. This method was also 
successfully used to tune the interatomic interaction be- 
tween atoms in the optical lattices, see, e.g., [a. [24 |48| . 

However, the correct theoretical treatment of a mag- 
netic Feshbach resonance requires the solution of a mul- 
tichannel problem that is numerically demanding due to 
the different length scales involved. Within an optical lat- 
tice the influence of the confining potential on the multi- 
channel problem has also to be properly considered (see 
|49j and references therein). To perform such a study 
within the present algorithm appears prohibitively dif- 
ficult, at least with the computer resources at our dis- 
posal. On the other hand, the variation of the interac- 
tion strength (characterized in the trap- free situation and 
at the zero-energy limit by a sc ) can be mimicked within 
single-channel models. In this case, some parameter in 
the rel. motion Hamiltonian is varied in such a fashion 
that resonant behavior occurs. Whenever the manipula- 
tion leads to a shift of a very weakly bound state into the 
dissociation continuum, resonant behavior (divergence of 
a sc ) is observed. Examples for such artificially obtained 
single-channel resonances include the variation of van der 
Waals coefficients [501, the reduced mass [23f | , th e inner- 
wall of the molecular interaction potential [26[, or the 
local modification of the Born-Oppenheimer curve at in- 
termediate distances [51J . A comparison of these different 



procedures and the full multi-channel treatment is pro- 



vided in [5l| . As is shown in [28[ a better and in fact 
almost perfect model for a multi-channel Feshbach reso- 
nance can be obtained with a two-channel model which 
appears more realistic with respect to a possible im- 
plementation within the present approach than the full 
multi-channel Hamiltonian. 



B. Basis-set considerations 

Due to the choice of spherical cm. and rel. coordi- 
nates all integrals could be reduced to products of one- 
dimensional integrals that can be solved very accurately 
and efficiently. In fact, the angular integrations are per- 
formed analytically. Moreover, the Gaussian quadra- 
ture provides exact results for all radial integrals except 
the ones of the interatomic interaction. However, even 
the latter ones can be calculated to high precision using 
Gaussian quadrature. The compactness of the B splines 
leads to sparseness of the Hamiltonian matrices, since 
only few integrals involve two non-zero B splines. The 
compactness and thus the resulting band structure of the 
one-particle Hamiltonian matrices (of cm. and rel. mo- 
tion) is controlled by the order fc rG i. and fc c .m. of the B 
splines. A higher order leads to a broader band struc- 
ture, but it offers also a higher flexibility of the basis 
functions, since it corresponds to a polynomial of higher 
order. Thus less B splines are needed for a comparable 
result, if a higher order is adopted. Usually, the orders 
of fc re i. (fcc.m.) of 8 or 9 turn out to be the best compro- 
mise with respect to basis-set size and sparseness, but 
also with respect to numerical stability in view of the fi- 
nite precision in which floating-point numbers are stored 
in the computer. 

The two other parameters defining the radial i?-splinc 
basis are the number of B splines and their knot se- 
quence. Clearly, the computational efforts (number of 
integrals that have to be calculated and size of the matri- 
ces that have to be diagonalized) depend crucially on the 
number of B splines. Since more B splines are required 
for describing highly oscillatory wavef unctions, it is most 
efficient to use non-uniform knot sequences in which the 
B-splinc density is higher in the highly oscillatory regions 
of (radial) space. 

In the context of ultracold collisions the energetically 
low-lying cm. orbitals ip(R) are the ones that usually 
are of main interest. Since they possess only a small 
number of nodes (none in the lowest state), the demands 
on the B-spline basis are not too high. For the simplest 
case of a single-well lattice potential Nr = 50 to Nr = 
100 B splines (or even less) were found to be sufficient 
to obtain convergence for the ground and lowest-lying 
excited states of the cm. motion |26|]. Also the lowest- 
lying states in a more structured trap potential like the 
triple-well potential considered in [28| were satisfactorily 
treated with 70 B splines. A linear knot sequence is 
usually adequate in this case. 

The description of the rel. orbitals ip(f) is more de- 
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manding, since one is usually not interested in the 
lowest-lying, deeply bound molecular states, but in the 
most weakly bound states or the low-lying dissociative 
ones. The Born-Oppenheimer potentials of alkali-metal 
atom dimers support often a large number of bound 
states [23]. The very long-ranged, weakly bound states 
consist, therefore, of a highly oscillatory inner part (cov- 
ering the molecular regime and providing the orthogonal- 
ity to all lower lying bound states) and a rather smooth 
long-range part. Hence, it is advantageous to distribute 
a majority of B splines in the molecular range of the po- 
tential while they are a sparser distributed in the outer 
part. The distribution of the B splines is given by the 
knot sequence {r-j}, i = 1, 2, . . . specifying a continuous 
chain of segments on which the B-spline functions are de- 
fined. The choice of a combined linear and geometrical 
knot sequence [521 ] for describing the short and long-range 
parts, respectively, has proven to be very efficient also in 
the present context [2{|. The linear distribution of the 
knot sequence is given by 



r i+fc« 



i = l,2 7 ... 7 Nl in - -k 



rel. 



(89) 



where AT' 1 "' and p m i n . are the number of B splines in 
the linear interval and the origin of the linear interval, 
respectively. Furthermore, s is the linear step size, i.e., 
the distance between neighboring knots given as 



Pin 



Ai 



M ir 



fcrel. + 1 



(90) 



Owing to the steep inner-wall of the molecular poten- 
tial the wave function <fi(r) vanishes well before r = 0. 
Therefore, p m i n . is usually taken non-zero in order to save 
on the number of B splines. The exact value of p m in. de- 
pends on the potential of the considered electronic state. 
Note, k r points must be placed at both ends of the box in 
accordance with the definition of B splines on the knot 
sequence (e.g., r\ = r% = ••• = rk rcl = Pmin.)- The linear 
step s is taken as the scale factor for the geometric pro- 
gression to ensure the smooth distribution of B splines at 
the border of the linear and geometric zones. In the geo- 
metrical knot sequence the separation of the knot points 
increases (according to a geometric series) with increas- 
ing distance. It is given by 



r- i+A riin. = piin. + sq 



i = l,...,N? 



(91) 



where A^ co - is the number of B splines in the geomet- 
ric interval and q is the common ratio for the geometric 
progression defined as 



1 



Pbo 



Aii 



N^ - 1 



(92) 



The last parameter is the maximum value of angular 
momentum l max used and thus the number of angular 
basis functions. Clearly, a more spherical-like lattice po- 
tential needs less angular momenta for representing the 



orbitals. The worst case is a highly anisotropic lattice 
geometry, since the spherical-harmonic basis converges 
extremely slowly in such a case. While a small value 
of ( max of 1 or 2 was sufficient to obtain converged or- 
bitals for an isotropic (cubic) single- well potential in [26[ 
in which among others also the experiment in |24( was 
successfully modeled and thus realistic parameters were 
adopted. On the other hand, the much more anisotropic 
triple- well potential considered in [28| required l max = 32 
for converged orbitals. 



C. Example of a very anisotropic trap 

Since a highly anisotropic lattice geometry provides 
a great challenge to our approach, it is important to 
demonstrate that even such a problem can be handled 
satisfactorily. Another motivation is the present interest 
in ultracold atomic systems of reduced dimensionality. 
Using strong confinement along one or two orthogonal 
directions, quasi one- or two-dimensional structures can 
nowadays be produced in the lab. Such systems show re- 
markable quantum properties not encountered in three 
dimensions. A major experimental breakthrough was 
the realization of To nk- Girardea u g ases of Bosons with 
strongly repulsive interaction [5j,[5J] . Being placed in ID 
and repelling each other Bosons are hindered from occu- 
pying the same position in space. This mimics the Pauli 
exclusion principle for Fermions, causing the Bosonic par- 
ticles to exhibit Fermionic properties. Another peculiar 
property of reduced dimensionality systems is the occur- 
rence of confinement-induced resonances that were re- 
cently experimentally observed [55l . I56J . Confinement- 
induced ID Feshbach resonances reachable by tuning 
the ID coupling constant via 3D Feshbach scattering 
resonances occur for both Bose gases [57| and spin- 
aligned Fermi gases [58|. Near a confinement-induced 
resonances, the effective ID interaction is very strong, 
leading to strong short-range correlations, breakdown of 
effective field theories, and emergence of a highly corre- 
lated ground state. Although confinement-induced reso- 
nances were originally predicted to occur already when 
only the rel. motion is considered [53, [59|, |60|, it was 
recently shown and furthermore verified using the here 
presented numerical approach that the ones observed in 
|55j are in fact caused by the coupling of the cm. and 
rel. degrees of freedom 61 1. 

A typical system to study effects of low dimension- 
ality is two atoms confined in a quasi ID cigar-shaped 
harmonic trap. We adopt this two-body setup to probe 
the quality of the present computational method in an 
extreme case of the strong anisotropic confinement. As 
was already mentioned, if the interparticle interaction 
between atoms is modeled by the (^-function pseudopo- 
tential the problem can be solved analytically. Therefore, 
the numerical results can be compared with the predic- 
tions of the pseudopotcntial model. One-dimensional ex- 
ternal confinement is prepared by setting a strong har- 
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monic frequency in cither of two spatial directions, c. g., 
lu x = ujy S> w 2 where u>i stands for the harmonic oscilla- 
tor frequency along the spatial direction i. It is worth- 
while to note that the analytical solution is only known 
for this special case of the single anisotropy [62(. This 
emphasizes the need for numerical approaches to cal- 
culate spectra in totally anisotropic confinement where 
u x ^ uj y > u g . 

The Hamiltonian of two identical particles in our har- 
monic trap separates in rel. and cm. motion. The cm. 
spectrum reduces to the one of a simple harmonic oscil- 
lator. This simplifies the problem to calculate the rel. 
spectrum. Hence, it is sufficient to solve numerically 
Eq. ([50)1 only. In the following two Bosonic 7 Li atoms 
arc considered that are placed in the prolate trap with 
uj x — uj y = 10 lo z and interacting via the potential of the 
a 3 £+ electronic triplet state. This electronic state has 
the advantage of supporting a small number of bound 
states. Therefore, a smaller number of B splines is re- 
quired to reproduce the ipi(r) functions. The numeri- 
cal data for the Born-Oppenheimer potential curve are 
adopted from [231 ]. In order to achieve converged results 
for the first few states Z max = 30, a box size of approx- 
imately pbox ~ 5 x 10 4 ao, where a$ is the Bohr radius, 
and N r = 100 B splines of the order fc rc i. = 8 were used. 
Half of the B splines (N l r in - = 50) are distributed lin- 
early according to Eq. (|89p over the small interatomic 
distance of pu n . = 15 clq to reproduce the highly oscilla- 
tory structure of v?i( r ) i n this region. Furthermore, for 
our calculations we can safely choose p m i n . = 2ao- The 
remaining B splines (N!° co - = 50) are distributed in an 
ascending geometric progression according to Eq. (|91[) 
over the residual interval. Finally, different values of the 
interaction strength are obtained using a single-channel 
approach by a smooth variation of the inner wall of the 
molecular interaction potential [5l[ . Figure 0] shows the 
calculated spectrum of the rel. Hamiltonian together with 
the analytical prediction of the pscudopotential approx- 
imation. The spectral curves are plotted as functions of 
dx/dsc, where d x = y/h/(fj,u} x ) is the harmonic oscillator 
length in transversal direction. As is evident from Fig. [4j 
the first four trap states and the bound state match per- 
fectly with the analytical prediction of the pscudopoten- 
tial model. However, high lying states show deviations. 

In order to explain the increasing mismatch found for 
higher lying states and to give an impression of the con- 
vergence behavior of the present approach, energy spec- 
tra for different values of the angular momentum arc 
calculated. For the calculations d x /a sc = 1.46 is cho- 
sen. Close to this value a confinement-induced resonance 
should occur [531 motivating this choice. The results of 
the calculations are shown in Fig. [5] As is seen from this 
figure, at Z max = 30 the slopes of the energy curves are ap- 
proximately zero, especially for the first four trap states. 
This indicates that convergence is achieved with respect 
to ! mal . Figure [5] clearly demonstrates that the high ly- 
ing states require larger values of the angular momentum 
for convergence. This can now as well be concluded from 




FIG. 4: Energy spectrum of two Bosonic atoms confined in a 
harmonic trap of anisotropy u) x /u) z = 10. The numerical cal- 
culation (blue dots) is compared to the analytical prediction 
of the pseudopotential approximation (red lines). 



Fig.m 




FIG. 5: The energy of the first five trap states for different 
values of l max in the strongly repulsive regime for d x /a sc = 
1.46. The deeper lying bound state is not shown here, because 
it is already converged for Z max = 5. (To guide the eye, the 
discrete points of the numerical calculation are connected by 
continuous lines.) 

While convergence can be extended to higher energies 
or higher anisotropies by increasing l max , the compu- 
tational efforts increase adding new angular momenta. 
However, as the precision is only limited by hardware ca- 
pacities, this example demonstrates the applicability of 
the approach even for extreme setups. 



D. Exact diagonalization (configuration 
interaction) 

The demands of the exact diagonalization also known 
as configuration interaction (CI) depend evidently again 
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on the physical system under consideration. While the 
atom-atom interaction is very efficiently handled even 
for strongly interacting atoms and basically contained in 
the rel. orbitals, the CI takes care of the cm. and rel. 
coupling. Thus its convergence depends on the strength 
of this coupling. Correspondingly, convergence is much 
faster for weakly coupled systems. If all configurations 
that can be built with the cm. and rel. orbitals are in- 
cluded, full CI is achieved. While being well defined, full 
CI scales very inconveniently with the number of orbitals. 
Therefore, it is in practice in most cases more advanta- 
geous to include only a limited number of the possible 
configurations. For example, the inclusion of the orbitals 
of very high energy does usually not lead to a notice- 
able improvement of the low lying states. Therefore, an 
energy cut-off may be introduced in the orbital selec- 
tion. If the interest is mainly on those states that are 
weakly bound or dissociative, but close to the dissocia- 
tion threshold, the deeply bound molecular states can be 
omitted from the CI configurations [26|]. Note, however, 
that a good representation of those states in the calcu- 
lation of the rel. orbitals is nevertheless important, since 
otherwise the rel. orbitals of the weakly bound states are 
not well described and this can in practice not be com- 
pensated by the CI calculation. In the worst case the 
calculation of the rel. orbitals provides too few bound 
states and then the nodal structure of the weakly bound 
states is evidently wrong. 



Although the coupled Hamiltonian matrices are usu- 
ally much larger than the uncoupled ones and are there- 
fore harder to diagonalize, there are two points easing 
the treatment of the coupled ones. First, the eigenvalue 
problem to be solved in the CI step is a standard one 
(meaning that the the basis is orthonormal) , while gen- 
eralized ones occur in the orbital calculations. Second, 
in many cases only a relatively small number of CI states 
is required. Therefore, it is possible to use Lanzcos or 
Davidson type diagonalization routines that (iteratively) 
provide a small number of eigenstates within a specific 
energy interval. Here we adopted the Davidson-based di- 
agonalization routine JADAMILU [63[ that is especially 
designed for the efficient diagonalization of large sparse 
matrices. 



Finally, it should be noted that the choice of expressing 
all wavefunctions in rel. and cm. coordinates is advanta- 
geous for computational reasons, but often not very help- 
ful in the interpretation of the results. Especially in the 
case of multi-well potentials the obtained wavefunctions 
and densities are often very complicated. Therefore, a 
further code was implemented that allows the application 
of the inverse coordinate transform from cm. -rel. coordi- 
nates to the absolute ones according to fi — R + [iir and 
?2 = R — /xiT to the wavefunctions, especially to ^(R, r) 
which allowed a much easier physical interpretation, c g., 
in El. 



VI. SUMMARY AND OUTLOOK 

An approach that allows the full numerical description 
of two ultracold atoms in a finite orthorhombic 3D opti- 
cal lattice is presented. The coupling between center-of- 
mass and relative motion coordinates is incorporated in 
a configuration-interaction manner and hence the full 6D 
problem is solved. An important feature is the use of real- 
istic interatomic interaction potentials adopting, e. g., nu- 
merically provided Born-Oppcnhcimer curves. The use 
of spherical harmonics together with B splines as basis 
functions and the expansion of the trap in terms of spher- 
ical harmonics leads to an analytical form of the matrix 
elements, except those of the interparticle interaction, 
if the interaction potential is defined only numerically. 
The sparseness of the Hamiltonian matrices due to the 
use of the compact radial B-spline basis is considered 
explicitly. This makes the method computationally effi- 
cient. Additionally, the lattice symmetry and a possible 
indistinguishability of the atoms (Bosonic or Fermionic 
statistics) is considered explicitly. This simplifies the cal- 
culations and helps to interpret the solutions. 

The here presented approach has already proven its 
applicability by considering the influence of the anhar- 
monicity in a single site of an optical lattice in [261 ] where 
a corresponding experiment |24| could be reproduced 
and analyzed. The validity range of the Bosc-Hubbard 
model together with an improved determination of the 
Bosc-Hubbard parameters was investigated considering 
a triple- well potential in [231 ■ 

The implemented approach was formulated in a rather 
general way in order to allow extensions in various di- 
rections. Since the optical-lattice potential is (via the 
Taylor expansion) expressed as a superposition of poly- 
nomials, it is rather straightforward to consider other 
than pure sin or cos 2 potentials, as long as they can be 
represented with (a reasonable number of) polynomials. 
This includes tilted lattices and supcrlattices. Care has, 
however, to be taken that the orthorhombic symmetry 
is either preserved, or new symmetry rules have to be 
implemented. 

Substitution of the numerical interatomic potential by, 
e.g., the Coulomb potential allows to describe either 
two electrons or an exciton in quantum-dot atoms or 
molecules. It is furthermore planned to implement non- 
isotropic interatomic potentials like dipole-dipole interac- 
tions as they are of interest for Cr or Rydberg atoms and 
for heteronuclear diatomic molecules. Finally, an exten- 
sion of the approach in the direction of time-dependent 
problems (with time dependent lattice or interatomic in- 
teraction potentials) is currently under way. 



Acknowledgments 

The authors thank Prof. Matthias Bollhofer for pro- 
viding the code JADAMILU and helpful comments for 
its implementation and use. The authors are grateful 



18 



to the Deutsche Forschungsgemeinschaft (within Sonder- 
forschungsbereich SFB 450), the Fonds der Chemischen 
Industrie, and the Humboldt Centre for Modern Optics 



for financial support. S. Grishkcvich gratefully acknowl- 
edges the integrating project of the European Commis- 
sion (AQUTE) for the finantial support. 



[1] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. 

Wieman, and E. A. Cornell, Science 269, 198 (1995). 
[2] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van 

Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, 

Phys. Rev. Lett. 75, 3969 (1995). 
[3] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and 

P. Zoller, Phys. Rev. Lett. 81, 3108 (1998). 
[4] M. Greiner, O. Mandel, T. Esslinger, T. Hansen, and 

I. Bloch, Nature 415, 39 (2002). 
[5] M. Kohl, H. Moritz, T. Stoferle, K. Giinter, and 

T. Esslinger, Phys. Rev. Lett. 94, 080403 (2005). 
[6] D. R. Meacher, Contemp. Phys. 39, 329 (1998). 
[7] I. Bloch, Nat. Phys. 1, 23 (2005). 
[8] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, 

A. Sen, and U. Sen, Adv. Phys. 56, 243 (2007). 
[9] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 

80, 885 (2008). 

[10] J. Weiner, V. S. Bagnato, S. Zilio, and P. S. Julienne, 
Rev. Mod. Phys. 71, 1 (1999). 

[11] T. Loftus, C. A. Regal, C. Ticknor, J. L. Bohn, and D. S. 
Jin, Phys. Rev. Lett. 88, 173201 (2002). 

[12] C. A. Regal, C. Ticknor, J. L. Bohn, and D. S. Jin, Na- 
ture 424, 47 (2003). 

[13] J. Hubbard, Proc.R.Soc. A 276, 238 (1963). 

[14] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, 
and M. Greiner, Nature 472, 307 (2011). 

[15] T. Kraemer, M. Mark, P. Waldburger, J. G Danzl, 
C. Chin, B. Engeser, A. D. Lange, K. Pilch, A. Jaakkola, 
H.-C. Nagerl, et al., Nature 440, 315 (2006). 

[16] M. Valiente, D. Petrosyan, and A. Saenz, Phys. Rev. A 

81, 011601 (2010). 

[17] M. Valiente, M. Krister, and A. Saenz, Europhys. Lett. 

92, 10001 (2010). 
[18] B. Chatterjee, I. Brouzos, S. Zollner, and P. Schmelcher, 

Phys. Rev. A 82, 043619 (2010). 
[19] D. DeMille, Phys. Rev. Lett. 88, 067901 (2002). 
[20] D. Jaksch and P. Zoller, Adv. Phys. 315, 52 (2004). 
[21] P.-I. Schneider and A. Saenz, Quantum computation 

with ultracold atoms in a driven optical lattice (2011), 

arXiv: 1103.4950. 
[22] E. L. Bolda, E. Tiesinga, and P. S. Julienne, Phys. Rev. A 

71, 033404 (2005). 
[23] S. Grishkevich and A. Saenz, Phys. Rev. A 76, 022704 

(2007). 
[24] C. Ospelkaus, S. Ospelkaus, L. Humbert, P. Ernst, 

K. Sengstock, and K. Bongs, Phys. Rev. Lett. 97, 120402 

(2006). 
[25] F. Deuretzbacher, K. Plassmeier, D. Pfannkuche, 

F. Werner, C. Ospelkaus, S. Ospelkaus, K. Sengstock, 

and K. Bongs, Phys. Rev. A 77, 032726 (2008). 
[26] S. Grishkevich and A. Saenz, Phys. Rev. A 80, 013403 

(2009). 
[27] J. R. Armstrong, N. T. Zinner, D. V. Fedorov, and A. S. 

Jensen, J. Phys. B 44, 055303 (2011). 
[28] P.-I. Schneider, S. Grishkevich, and A. Saenz, 

Phys. Rev. A 80, 013404 (2009). 



[29] G J. Milburn, J. Corney, E. M. Wright, and D. F. Walls, 

Phys. Rev. A 55, 4318 (1997). 
[30] M. Anderlini, J. Sebby-Strabley, J. Kruse, J. V. Porto, 

and W. D. Phillips, J. Phys. B 39, S199 (2006). 
[31] S. Foiling, S. Trotzky, P. Cheinet, M. Feld, R. Saers, 

A. Widera, T. Midler, and I. Bloch, Nature 448, 1029 

(2007). 
[32] P. Cheinet, S. Trotzky, M. Feld, U. Schnorrberger, 

M. Moreno-Cardoner, S. Foiling, and I. Bloch, 

Phys. Rev. Lett. 101, 090404 (2008). 
[33] J. Sebby-Strabley, M. Anderlini, P. S. Jessen, and J. V. 

Porto, Phys. Rev. A 73, 033605 (2006). 
[34] S. Trotzky, P. Cheinet, S. Foiling, M. Feld, U. Schnor- 
rberger, A. M. Rey, A. Polkovnikov, E. A. Demler, 

M. Lukin, and I. Bloch, Science 319, 295 (2008). 
[35] J. A. Stickney, D. Z. Anderson, and A. A. Zozulya, 

Phys. Rev. A 75, 013608 (2007). 
[36] D. Blume and C. H. Greene, Phys. Rev. A 65, 043613 

(2002). 
[37] T. Busch, B.-G Englert, K. Rzazewski, and M. Wilkens, 

Found. Phys. 28, 549 (1998). 
[38] Z. Idziaszek and T. Calarco, Phys. Rev. A 71, 050701(R) 

(2005). 
[39] W. Kohn, Phys. Rev. 115, 809 (1959). 
[40] I. Gradshteyn and I. Ryzhik, Table of Integrals, Series, 

and Products (Academic Press, 2007). 
[41] S. Zollner, H.-D. Meyer, and P. Schmelcher, Phys. Rev. A 

78, 013621 (2008). 
[42] C. de Boor, A Practical Guide to Splines (Springer, New 

York, 1978). 
[43] H. Bachau, E. Cormier, P. Decleva, J. E. Hansen, and 

F. Martin, Rep. Prog. Phys. 64, 1815 (2001). 
[44] J. Rasch and A. C H. Yu, SIAM J. Sci. Comput. 25, 1416 

(2003). 
[45] D. Pinchon and P. E. Hoggan, Int. J. Quant. Chem. 107, 

2186 (2007). 
[46] F. C. von der Lage and H. A. Bethe, Phys. Rev. 71, 612 

(1947). 
[47] C. A. Regal and D. S. Jin, Phys. Rev. Lett. 90, 230404 

(2003). 
[48] T. Stoferle, H. Moritz, K. Giinter, M. Kohl, and 

T. Esslinger, Phys. Rev. Lett. 96, 030401 (2006). 
[49] P.-I. Schneider, Y. V. Vanne, and A. Saenz, Phys. Rev. A 

83, 030701 (2011). 
[50] E. Ribeiro, A. Zanelatto, and R. Napolitano, 

Chem. Phys. Lett. 390, 89 (2004). 
[51] S. Grishkevich, P.-I. Schneider, Y. V. Vanne, and 

A. Saenz, Phys. Rev. A 81, 022719 (2010). 
[52] Y. V. Vanne and A. Saenz, J. Phys. B 37, 4101 (2004). 
[53] M. Girardeau, J. Math. Phys. 1, 516 (1960). 
[54] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Foiling, 

I. Cirac, G V. Shlyapnikov, T. W. Hansch, and I. Bloch, 

Nature 429, 277 (2004). 
[55] E. Haller, M. J. Mark, R. Hart, J. G Danzl, L. Reich- 

sollner, V. Melezhik, P. Schmelcher, and H.-C. Nagerl, 

Phys. Rev. Lett. 104, 153203 (2010). 



19 



[56] B. Frohlich, M. Feld, E. Vogt, M. Koschorreck, W. Zw- 

erger, and M. Kohl, Phys. Rev. Lett. 106, 105301 (2011). [61] 

[57] M. Olshanii, Phys. Rev. Lett. 81, 938 (1998). 

[58] B. E. Granger and D. Blume, Phys. Rev. Lett. 92, 133202 [62] 
(2004). [63] 

[59] D. S. Petrov, M. Holzmann, and G. V. Shlyapnikov, 
Phys. Rev. Lett. 84, 2551 (2000). 

[60] T. Bergeman, M. G. Moore, and M. Olshanii, 



Phys. Rev. Lett. 91, 163201 (2003). 

S. Sala, P.-I. Scheider, and A. Saenz, Confinement- 
induced resonances revisited (2011), arXiv:1104.1561. 
J.-J. Liang and C. Zhang, Phys. Scr. 77, 025302 (2008). 
M. Bollhofer and Y. Notay, Comp. Phys. Comm. 177, 951 
(2007). 



